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ABSTRACT 



^ : 



£f*} . Aims. We present the results of hydrodynamic simulations of the formation and subsequent orbital evolution of giant planets embedded in a 
' circumbinary disc. The aim is to examine whether or not giant planets can be found to orbit stably in close binary systems. 
Methods. We performed numerical simulations using a grid-based hydrodynamics code. We assume that a 20 M© core has migrated to the 
edge of the inner cavity formed by the binary where it remains trapped by corotation torques. This core is then allowed to accrete gas from 
the disc, and we study its orbital evolution as it grows in mass. For each of the two accretion time scales we considered, we performed three 
simulations. In two of the three simulations, we stop the accretion onto the planet once its mass becomes characteristic of that of Saturn or 
O ' Jupiter. In the remaining case, the planet can accrete disc material freely in such a way that its mass becomes higher than Jupiter's. 
|M Results. The simulations show different outcomes depending on the final mass m p of the giant. For m p = 1 M s (where M s is Saturn's mass), 
C/} . we find that the planet migrates inward through its interaction with the disc until its eccentricity becomes high enough to induce a torque 
, ^ , reversal. The planet then migrates outward, and the system remains stable on long time scales. For m p > IMj (where Mj is Jupiter's mass) 
we observed two different outcomes. In each case the planet enters the 4:1 resonance with the binary, and resonant interaction drives up the 
H ■ eccentricity of the planet until it undergoes a close encounter with the secondary star, leading to scattering. The result can either be ejection 
^ ' from the system or scattering out into the disc followed by a prolonged period of outward migration. These results suggest that circumbinary 
planets are more likely to be quite common in the Saturn-mass range. Jupiter-mass circumbinary planets are likely to be less common because 
■ of their less stable evolution, but if present are likely to orbit at large distances from the central binary. 
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Ctf ■ 30 reside in binary or multiple-star systems, with most of circumbinary planets have not yet been observed in binary 

them orbiting one stellar component in so-called S-type orbits systems composed of two main sequence stars . 

(Eggenbergeretal. 2004; Mugraueretal. 2007). The majority Seyeral circumbinary discSi however , have been detected 

of binary stars hosting planets have orbital separation a b > arQund spectroscopic binaries like DQ TaU; AK S co and GW 

100 AU. However, there are a few cases of short period binary Qri In QQ ^ the circumbinary disc has been directly imaged 

systems with a b ~ 20 AU like Ghese 86, y Cephei and HD and has reyealed the presence of an inner disc cayity due t0 thg 

41004 in which planets have been discovered to orbit at 1-2 Mgi torques exerted by the cental binary (Dutrey et al 1994) 

AU from the primary (Eggenberger et al. 2004, Mugrauer & The existence of these circumbinary discs combined with the 

Neuhauser 2005). As suggested by studies of the long-term fact ^ „ 5Q % of solar . type stars are members of binaries 

stability of planets in binary systems (e.g. Holman & Wiegert (Duquennoy & Mayor 1991) suggests that circumbinary 

1999), close binaries with a b ~ 1 AU can, in principle, harbour planets may be CQmmon provided that planet formation can 

planets evolving on a P-type orbit which encircles the two occur inside such discs, 
components of the binary system. One such circumbinary 
planet with mass of m p = 2.5 Mj has been detected orbiting 

at 23 AU from the radio pulsar binary PSR 1620-26. Another To date, several theoretical studies focused on planet 

formation in close binary systems indicate that planetesimal 

Send offprint requests to: A. Pierens accretion should be possible within circumbinary discs. 



1. Introduction with mass of m p = 2.44 My was found evolving around a 

system which consists of the star HD 202206 and its 17.4 Mj 
Among the approximately 260 extrasolar planets known at brQwn dwarf companion (Udry et aL 20 02). Because short 
the time of writing (e.g. http://exoplanet.eu/), about period binaries are often rejected from observational surveys , 
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Moriwaki & Nakagawa (2004) found that planetesimals can 
grow in gas-free circumbinary discs only in regions farther out 
that ~ 13 AU from a binary with orbital separation = 1 AU, 
eccentricity = 0.1 and mass ratio = 0.2. The influence 
of gas drag was studied recently by Scholl et al. (2007) who 
showed that for a binary with the same parameters, the eccen- 
tricity damping provided by the disc can enable planetesimal 
accretion to occur in regions located below ~ 4 AU from the 
central binary. The later stages of planet formation in which 
Earth-mass planets form by accumulation of embryos was 
investigated by Quintana & Lissauer (2006). These authors 
found that planetary systems similar to those around single 
stars can be formed around binaries, provided that the ratio 
of the binary apocentre distance to planetary orbit is < 0.2. 
In general binaries with larger maximum separations lead to 
planetary systems with fewer planets. 

Recently, the evolution of Earth-mass bodies embedded in a 
circumbinary disc and undergoing type I migration because 
of disc torques (e.g. Ward 1997) was examined by Pierens & 
Nelson (2007) (hereafter referred to as Paper I). In this work, it 
was found that the inward drift of a protoplanet can be stopped 
near the edge of the cavity formed by the binary. Such an effect 
arises because in this region, the gradient of the disc surface 
density is such that the planet experiences strong positive 
corotation torques which can eventually counterbalance the 
negative differential Lindblad torque, thereby leading to the 
halting of migration (Masset et al. 2006). In a subsequent 
paper, Pierens & Nelson (2008a) extended this work by inves- 
tigating the issue of how multiple protoplanets interact with 
each other if they form at large distance from the binary and 
successively migrate toward the cavity edge. The simulations 
performed by Pierens & Nelson (2008a) of pairs of planets 
interacting with each other indicated different outcomes such 
as resonant trapping or orbital exchange, depending on the 
ratio between the masses of the planets. Interestingly, in 
simulations involving more than two planets, planetary growth 
resulting from scattering and collisions between protoplanets 
was found to be feasible. This implies that giant cores might 
be formed in circumbinary discs, resulting eventually in a gas 
giant planet orbiting near the cavity edge. 
The orbital evolution of Jovian mass planets embedded in 
circumbinary discs was studied by Nelson (2003). This work 
showed that giant planets undergoing type II migration (e.g. 
Lin & Papaloizou 1993; Nelson 2000) are likely to enter the 
4:1 resonance with the binary. In that case, the subsequent 
evolution of the giant can be twofold, depending on whether 
or not the resonance is stable. For systems in which the 4:1 
resonance is stable, the planet remains near or at the resonance. 
However, it appears that there is a finite probability for the 
system to be unstable, in which case the planet can be ejected 
from the system due to close encounters with the binary. 

In this paper, we extend the work of Nelson (2003) by 
studying the whole evolution of a circumbinary planet during 
its growth from a core into a gas giant. To address this issue, 
we consider a scenario in which a 20 M core initially trapped 
at the edge of the cavity can slowly accrete gas from the 
disc. As the planet grows, the onset of non-linear effects 



as well as gap formation can exclude gas material from the 
coorbital region, thereby cancelling the effects of corotation 
torques. This can subsequently lead to the planet migrating 
inward again, until it becomes trapped eventually in a mean 
motion resonance with the binary. Here, we present the results 
of hydrodynamic calculations aimed at simulating such a 
scenario. In particular, we want to examine how the final 
outcome of the system depends on the accretion rate onto the 
planet as well as on the final mass of the giant. Interestingly, 
the results of the simulations indicate that only Saturn-mass 
giant planets can evolve stably in circumbinary discs. Most of 
the calculations of embedded giants with masses of m p > 1 Mj 
resulted in close encounters between the planet and the binary, 
leading eventually to the planet being completely ejected from 
the system. 

This paper is organized as follows. In Section 2, we de- 
scribe the hydrodynamical model. The results of the simula- 
tions are discussed in Section 3. We finally summarise and 
present our conclusions in Section 4. 

2. Hydrodynamical model 

2.1. Numerical method 

We consider a 2D disc model in which all physical quanti- 
ties are vertically averaged and we work in polar coordinates 
(r, 0) with the origin located at the centre of mass of the bi- 
nary. The equations governing the disc evolution as well as the 
equations of motion for the binary plus planet system can be 
found in Paper I. The equations for the disc are solved using the 
hydrocode Genesis, which has been tested extensively against 
other codes (De Val-Borro et al. 2006), and employs a second 
order numerical scheme based on the monotonic transport al- 
gorithm (Van Leer 1977). Included in this code is a fifth-order 
Runge-Kutta integrator (Press et al. 1992) used to compute the 
evolution of the planet and binary orbits. 
As in Paper I, we use N r = 256 radial grid cells uniformly dis- 
tributed between r in = 0.5 and r out = 6 and = 380 azimuthal 
grid cells. We adopt also the same computational units in which 
the mass of the binary is M+ = 1 , the gravitational constant is 
G = 1 and the radius r = 2 in the computational domain corre- 
sponds to 5 AU. The unit of time is D" 1 = JGM*/a 3 b , where 
ab = 0.4 is the initial binary separation. In the following, we 
report our results in units of the initial orbital period of the bi- 
nary P = ItiQT 1 . 

In this paper, the planet is allowed to accrete gas from the disc. 
Accretion by the protoplanet can be modelled by removing at 
each time-step a fraction of the gas located inside the Roche 
lobe of the planet and then adding the corresponding amount 
of matter to the mass of the planet (e.g. Kley 1999; Nelson et 
al. 2000). Here, the removal rate is chosen such that the accre- 
tion time scale onto the planet is t acc = td yn , where td yn is the 
orbital period of the planet. This corresponds to the maximum 
rate at which the planet can accrete gas material (Kley 1999). In 
order to examine how the final outcome of the system depends 
on the accretion rate, we have also performed simulations with 

tacc = 10 tdy n . 
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Fig. 1. This figure shows the initial disc surface density profile. 
The initial position of the planet with mass of m p = 20 M© is 
represented by a black circle. 



When calculating the gravitational force between the disc 
and planet we adopt a gravitational softening parameter b = 
0.6//, where H is the disc height at the planet location. Material 
contained within the planet Hill sphere does not contribute to 
the gravitational acceleration experienced by the planet. 

2.2. Initial conditions 

In Paper I, we presented the results of simulations of proto- 
planets with masses of m p = 5, 10 and 20 M© embedded in 
circumbinary discs. We found that in each case, the inward 
migration of the protoplanet is halted due to the action of 
corotation torques which operate strongly at the cavity edge. 
Here, we extend the model presented in this earlier paper and 
examine how such a trapped protoplanet evolves as it accretes 
gas from the disc and grows to become a giant planet. In order 
to investigate this issue, we restarted one of the simulations 
presented in paper I for which m p = 20 M© at a point in time 
when the planet is trapped at the cavity edge. Therefore, in this 
new series of simulations, the protoplanet evolves initially on 
an orbit with a p ~ 1.2 and e p ~ 0.02, which are the values for 
the semi-major axis and eccentricity that a 20 M© body finally 
attains once its migration has been stopped (see Paper I). 
The initial semi-major axis and eccentricity of the binary are 
at ~ 0.39 and ~ 0.08 respectively. These values correspond 
to the ones that an initially circular binary with - 0.4 
and mass ratio = 0.1 ultimately attains as it interacts 
with a circumbinary disc. In Paper I, we indeed showed that 
the evolution outcome of such a system is an equilibrium 
configuration for which the disc structure as well as the binary 
eccentricity remain unchanged. Interestingly, from the time 
this quasi-steady state is reached, we found that the apsidal 
lines of the disc and binary are almost aligned. 

The disc model used in this work is the same as that of 
Paper I. Accordingly, the aspect ratio is constant and equal to 
H/r = 0.05. The initial disc surface density profile is presented 
in Fig. [U For numerical reasons (see Section 12.31) , we use a 
low-density region from r = 4 to r - 6 . From r ~ 2 to r - 4, 



the surface density is S(r) = r" 1/2 where So is chosen such 
that the unperturbed disc would contain 0.01 M© inside 10 AU 
(we assume that the initial binary separation corresponds to 1 
AU). We further note that such a density profile corresponds 
to the equilibrium density for an accretion disc with a constant 
kinematic viscosity (e.g. Gunther & Kley 2002). In the inner 
parts, the disc presents an inner cavity which is created self- 
consistently from simulations of binary-disc interactions (see 
Paper I). It is worthwhile to notice that here, the onset of non- 
linear effects due to the presence of a 20 M© initially located at 
r ~ 1.2 can modify slightly the gap profile. 
We model the disc turbulent viscosity using the "alpha" pre- 
scription for the effective kinematic viscosity v = ac s H 
(Shakura & Sunyaev 1973). In Paper I, we set a = 10" 4 be- 
cause using larger values caused the binary separation to de- 
crease too rapidly to allow an equilibrium configuration to be 
obtained (see Paper I for details). In this work however, we 
examine the evolution of giant planets undergoing type II mi- 
gration. Because in that case the migration rate of the planet 
is controlled by the disc viscous evolution, we decided here to 
use a more realistic a value, which probably lies in the range 
10" 3 - 10" 2 in real circumstellar dies. In order to obtain a disc 
model in which a = 10" 3 , the calculations were started with a 
increasing slightly from 10" 4 to the desired value in ~ 1500 bi- 
nary orbits. Although not sufficient for a new equilibrium con- 
figuration to be established, this gives a sufficient time to the 
binary plus disc system to readjust to the new a value. 

2.3. Boundary conditions 

In order to avoid any wave reflection at the outer edge of the 
computational domain, we impose a low-density region be- 
tween r = 4 and r = 6 using a taper function. 
At the inner boundary, we model the accretion onto the central 
star by setting the radial velocity in the inner ghost zones to 
v r = fiVriXin), where v r (r ;w ) = -3v/2r;„ is the gas drift velocity 
at the disc inner edge due to viscous evolution and where ft is 
a free parameter. Following Pierens & Nelson (2008b), we set 
J3 = 5 in the simulations presented here. 

3. Results 

For each value of the accretion rate onto the planet we consider, 
we have performed three simulations which differ in the final 
mass attained by the giant. In the first and second calculations, 
we prevent further growth of the planet once its mass becomes 
m p = 1 M s (where M s is Saturn's mass) and m p = 1 Mj 
respectively. In the last run however, the final mass of the planet 
is not restricted and we allow the latter to accrete gas freely 
from the disc. The mass of the planet as a function of time is 
presented in Fig. [2] for each run. 

3.1. Evolution of Saturn-mass planets 

For both values of t acc , the semi-major axis evolution of a 20 
M© body which grows to become a Saturn-mass planet is dis- 
played in the upper panel of Fig. [3] In the model with t acc = td yn , 
the planet mass reaches m p = 1 M s at t ~ 300 P, and this mass 
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Fig. 3. This figure shows the evolution of the system for models in which m p = I Ms. Here, the left and right panels correspond 
to t acc = tdyn and t acc = 10 td yn respectively. Top: evolution of the semi-major axis of the planet. Middle: evolution of the 
eccentricities for the planet (black), binary (red) and disc (green). Bottom: evolution of the longitudes of pericentre for the planet 
(black), binary (red) and disc (green). 
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Fig. 2. This figure shows the mass of the planet as a function 
of time for simulations in which t acc = td yn {upper panel) and 
tacc = 10 tdyn (lower panel). The dotted (resp. dashed) line 
corresponds to simulations in which the final mass of the planet 
is 1 M s (resp. 1 Mj). In simulations represented by the solid 
line, the final mass of the planet is not imposed. 



is attained at t ~ 2800 P in the model with t acc = 10 td yn (see 
Fig. [2]). The orbital evolution of the planet is, however, very 



similar in both cases and typically proceeds as follows. At the 
beginning of the simulation, gap formation due to the growth of 
the protoplanet can exclude material from the coorbital region, 
which results in the (negative) differential Lindblad torques 
being no longer couterbalanced by the (positive) corotation 
torques. The planet thus migrates inward again. As it migrates, 
the interaction with the central binary makes e p increase slowly, 
as observed in the middle panel of Fig. [3] which shows the evo- 
lution of the eccentricities of the binary, disc and planet for both 
values of t acc . We define the disc eccentricity ed by: 



where r max is set to r max = 3 and where e c is the eccentricity of 
the disc computed at the centre of each grid cell. This orbital 
element can be deduced by computing the eccentricity vector 
e = (v x h)/M* - r/r, where v and h are respectively the 
velocity and angular momentum of the disc at the cell centre. 

Once the planet eccentricity reaches e p ~ 0.05, which 
occurs at t ~ 500 P, the migration of the planet stops and 
then reverses. Previous work of Papaloizou & Larwood 
(2000) indicated that an eccentric low-mass protoplanet can 
experience net positive disc torques when e p ~ 1.1 H/R, owing 
to the latter crossing resonances in the disc that do not overlap 
the orbit at low-eccentricities. As pointed out by these authors, 



A. Pierens and R.P Nelson: On the formation and migration of giant planets in circumbinary discs 



4x1 G -5 I I I I 1.10 




-4x1 0" 5 I I I I I 0.95 

2.88000x10 4 2.88075x10 4 2.88150x10 4 2.88225x10 4 2.88300x10 4 

Time (binary orbits) 

Fig. 4. This figure shows, over a few planetary orbital periods, 
the evolution of the torques exerted by the disc on the planet 
(solid line) as well as the time evoution of the orbital the posi- 
tion of the planet r p (dashed line). 



this occurs because a planet evolving on a high eccentric orbit 
rotates more slowly than the disc at apocentre, which results 
in the outer disc exerting a positive torque on the planet. 
Although here the planet mass is in the Saturnian mass range, 
it appears that the observed migration reversal is caused by 
a similar phenomenon. In order to demonstrate that such an 
effect is at work here, we have plotted in Fig. [4] the evolution 
of both the disc torques and planet orbital position r p over a 
few orbital periods of the planet. In agreement with Papaloizou 
& Larwood (2000), the disc torques exerted on the planet are 
positive (negative) when the latter is at apocentre (pericentre). 
Clearly, there is a slight imbalance between the positive and 
negative torques, due to the planet orbiting at a cavity edge, 
which combined with the fact that the planet spends more time 
at apocentre, favours a net positive torque on the planet and 
outward migration. We note that in Fig. |H there is a slight 
phase shift between the two curves which corresponds to the 
time needed for the planet to create an inner (outer) wake at 
apocentre (pericentre). 

Whereas this reversed migration proceeds very slowly until 
t ~ 4.5 x 10 4 P, subsequent evolution involves exponential 
growth of the planet migration rate, which is characteristic 
of an episode of runaway migration (Masset & Papaloizou 
2003). Saturn-mass planets are known to be good candidates 
for such a migration regime in massive protoplanetary discs, 
due to their ability to create a coorbital mass deficit larger 
than their own mass. As noticed by Masset & Papaloizou 
(2003), outward runaway migration can eventually occur in 
discs with shallow surface density profiles, provided that the 
planet initially migrates with a significantly positive drift 
rate. Here, two effects contribute to make such a phenomenon 
possible. First, prior to the period of runaway migration, the 
planet migrates outward over ~ 2 x 10 4 binary orbits, which is 
clearly larger than the planet libration timescale. Second, the 
planet evolves in a region of strong positive surface density 
gradient such that the coorbital mass deficit increases as the 
planet migrates outward. This can be seen for example in Fig. 
[5] which shows, for the model in which t acc = td yn , a series of 
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surface density plots at different times. 

Interestingly, this period of outward runaway migration breaks 
at t ~ 8 x 10 4 P, which corresponds to the time when the planet 
passes through the 5:1 resonance with the binary. This can 
excite the planet eccentricity to such values that during the 
course of an orbit, the radial excursion of the planet can exceed 
its coorbital width, leading to the loss of the coorbital mass 
deficit. The fourth panel in Fig. \5\ displays the disc surface 
density just after the 5:1 resonance crossing. Comparing this 
panel with the second one which corresponds to a time when 
the planet is about to undergo runaway migration, it is clear 
that the coorbital region is no longer depleted, which can 
prevent runaway migration being sustained. This result is in 
broad agreement with the calculations performed by Masset & 
Papaloizou (2003) which indicate a similar tendency for the 
runaway migration to not be maintained. 

After the 5:1 resonance crossing, continuation of the runs 
indicates that the planets undergo slow outward migration for 
~ 5 x 10 4 binary orbits and then migrate inward again. Due to 
the very long run times required by these simulations, the final 
fate of the planets remains uncertain. Nevertheless, it seems 
likely that continued inward migration will proceed until the 
eccentricities of the planets are high enough for the disc torques 
to become positive, resulting eventually in the planets migrat- 
ing outward again. Thus we expect the long term evolution to 
consist of periods of inward followed by outward migration 
until disc dispersal, resulting in a stable Saturn-mass circumbi- 
nary planet. 

The evolution of the longitudes of pericentre of the binary, 
disc and planet is depicted in the lower panel of Fig. [51 As in 
Paper I, we compute the disc longitude of pericentre according 
to the following definition: 

Wd = r2n rr ; (2) 

where w c is the longitude of pericentre of the disc at the cen- 
tre of each grid cell. This orbital element is computed from 
the knowledge of the eccentricity vector using the relation 
cos(w c ) = e x /e. We note that Wd may be significantly different 
if the gravity of the disc was taken into account in the simula- 
tions. Indeed, it is well-known that one of the main effects of 
self-gravity is to induce prograde precession of the disc (e.g. 
Papaloizou 2002). The figure shows that the apsidal line of the 
planet is almost perfectly aligned with the ones of the disc and 
binary. Here, it is worthwhile noticing that despite the fact that 
the disc and binary precess at the same rate, the calculations 
show continued growth of the binary eccentricity due to reso- 
nant interactions between the disc and binary (see middle panel 
of Fig. [3]). This indicates that the binary plus disc system has 
not yet reached an equilibrium configuration where the eccen- 
tricity forcing has saturated. 

3.2. Evolution of Jupiter-mass planets 

For the two accretion time scales we consider, the orbital 
evolution of a planet with final mass m p = 1 Mj is illustrated 
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Fig. 5. This figure shows, for the model in which m p = 1 Ms and t acc = td yn , snapshots of the disc surface density at times shown 
above the plots. Here, the planet is represented by a white circle. 



in the upper panel of Fig. [6l In the model with t acc - td yn the 
planet reaches one Jupiter mass in ~ 5.5 x 10 3 binary orbits 
while it reaches this mass in ~ 1.4 x 10 4 orbits in the model 
with t acc = 10 tdyn (see Fig. O. The early evolution of the 
planet is similar to that described for Saturn-mass planets, 
involving inward migration of the planet and continued growth 
of its eccentricity. However, contrary to models in which 
m p - 1 M$, simulations with m p = 1 Mj resulted in the 
temporary formation of the 4:1 resonance between the planet 
and binary (we note that migration reversal occurs for the 
planets with m p = 1 Ms before they reach the 4:1 resonance). 



For the calculation with t a , 



tdyn, the time evolution of 



the resonant angles iff\ 9 fa, 1J/3 and ^4 associated with the 4:1 
resonance is displayed in Fig. [71 These are given by: 



(3) 



where (A p ) and cjjj (oj p ) are respectively the mean longitude 
and longitude of pericentre of the binary (planet). Once the 
resonance is established, which occurs at t ~ 3 x 10 3 P for 
this model, the resonant angle if/2 librates with low amplitude, 
thereby indicating that the planet is strongly locked into the 
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resonance. The evolution of the eccentricities of the planet, 
disc and binary is illustrated in the middle panel of Fig. 
[6) As expected, the resonant interaction with the binary is 
accompanied by a significant growth of the planet eccentricity 
e p , which increases up to e p ~ 0.15. At t ~ 8.5 x 10 3 P, this 
eccentricity growth causes the protoplanet to undergo a close 
encounter with the binary, leading subsequently to the planet 
being extracted from the resonance and scattered further out in 
the disc on a high-eccentric orbit with a p ~ 2. 

A similar outcome is observed in the simulation with 
tacc = 10 tdyn- Here the planet is scattered at t ~ 4 x 10 4 P 
and evolves subsequently on an orbit with significantly larger 
semimajor axis (a p ~ 1.6) and eccentricity (e p ~ 0.3). Since 
a close encounter is generally a chaotic dynamical process, 
these values are quite different from the ones observed in the 
simulation with t acc = td yn - In fact, these strongly depend on 
the encounter geometry and therefore may significantly differ 
from one simulation to the other. 

Prior to this close encounter however, the evolution of the 
planet differed noticeably from that corresponding to the 
model with t acc = td yn - It appears that in the calculation with 
tacc = 10 tdyn, the 4:1 resonance becomes rapidly undefined 
after its formation, which arises at t ~ 4 x 10 4 . This is because, 
compared with the run in which t acc - td yn , the planet mass 



A. Pierens and R.P Nelson: On the formation and migration of giant planets in circumbinary discs 



7 



t QCC =U 




1.0x10 3 
Time (binary orbits) 



t acc =10 t dyn 









£ .. nlllll(lpll J 

















10 4 8.0x1 
Time (binary 



3 1.0x10 s 

orbits) 




5.0x10 4 1.0x1 5 1.5x10 5 2.0x10 4 4.0x1 4 6.0x10 4 8.0x10 4 1.0x10 5 1.2x10 5 1.4x10 5 

Time (binary orbits) Time (binary orbits) 



Fig. 6. This figure shows the evolution of the system for models in which m p = 1 Mj. Here, the left and right panels correspond 
to t acc = tdyn and t acc = 10 td yn respectively. Top: evolution of the semi-major axis of the planet. Middle: evolution of the 
eccentricities for the planet (black), binary (red) and disc (green). Bottom: evolution of the longitudes of pericentre for the planet 
(black), binary (red) and disc (green). 



is significantly lower when the 4:1 resonance is established, 
leading to a weaker resonant locking. Once the resonance is 
broken, the planet is located just outside of the 4:1 resonance 
and migrates slightly outward until t ~ 1.7 x 10 4 P. Outward 
migration is induced by positive disc torques due to the 
planet eccentricity having reached e p ~ 0.15 after resonance 
breaking. Then, the planet migrates inward and approaches the 
4:1 resonance again, but its eccentricity remains sufficiently 
high for there to be a close encounter between the planet and 
central binary system leading to the planet being scattered by 
the binary. The planet and binary undergo multiple encounters 
during this phase, with significant changes in their orbital 
elements occuring (see middle-right panel of Fig. 0). The 
planet ends up orbiting interior to the 4:1 resonance where it 
migrates inward toward the 3:1 resonance. 

In both simulations, the evolution of the planet subsequent 
to this initial scattering is quite similar. In the following, we 
use the results of the simulation with t acc = td yn to discuss in 
more details how the evolution of the system proceeds after 
this event. At t ~ 10 4 P, interaction with the disc and central 
binary finally results in the planet settling into an orbit further 
out in the disc with eccentricity e p ^ 0.1-0.15. The planet 
now forms a gap in the disc and begins to migrate inward 
slowly under type II migration. The evolution of the disc and 
planet plus binary system for this run is presented in Fig. [8] 



which shows snapshots of the disc surface density at different 
times. The first panel corresponds to a time shortly before 
the formation of the 4:1 resonance while the three other ones 
display the state of the system just after the initial scattering. 
These show clearly that the inner disc is progressively lost 
through the inner boundary as a result of viscous evolu- 
tion, thereby leading to a reduction of the (positive) inner disc 
torques and consequently to the inward migration of the planet. 

At the end of the run with t acc = td yn , the final fate of the 
planet is still uncertain despite the very long time scale covered 
by the simulation. However, it is likely that the continued in- 
ward migration of the planet, combined with the growth of the 
eccentricities of both the binary and planet (see middle panel 
of Fig. [6]), will result in further close encounters and eventually 
in the planet being completely ejected from the system. Such 
an outcome is observed in the calculation with t acc = 10 td yn - 
In that case, this arises because inward migration of the planet 
causes the temporary formation of the 3:1 resonance with the 
binary at t ~ 1.35 x 10 5 P, which makes the planet eccentric- 
ity increase up to e p ~ 0.15. This leads to close encounters 
between the planet and the secondary star, resulting in the scat- 
tering and ejection of the planet from the system. 
The lower panel of Fig. [6] shows the evolution of the disc, 
binary and planet longitudes of pericentre for both models. 
Interestingly, we see that here the planet is aligned with nei- 
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Fig. 8. This figure shows, for the model in which m p = 1 Mj and t acc = tdyn, snapshots of the disc surface density at times shown 
above the plots. Here, the planet is represented by a white circle. 



ther the disc nor the binary. Therefore, relative to simulations 
with m p = I Ms, the probability of close encounters between 
the planet and binary is increased, thereby leading to a system 
which is more likely to become unstable. 
Lastly, we notice that the results of these runs are broadly con- 
sistent with previous hydrodynamical simulations of jupiter- 
mass planets embedded in circumbinary discs (Nelson 2003) 
which indicated that indeed, trapping into 4:1 resonance fol- 
lowed by a scattering through a close encounter with the binary 
is a possible outcome of such systems. 



3.3. Evolution of giants with m p > 1 Mj 

The evolution of accreting bodies with final masses m p > 1 Mj 
is depicted, for both values of the accretion time scale, in the 
upper panel of Fig. [9] With respect to the simulations presented 
in Sect. 13.21 here the calculations differ in that the planet can 
continue to accrete gas once its mass has attained m p = 1 Mj. 
Although subsequent evolution may differ, this implies that 
while m p < 1 Mj the evolution of the planet is the same as the 
planets considered in Sect. 13.21 Therefore, when discussing 
the results below, we consider only the evolution of the system 
from the time when the planet mass has reached m p = 1 Mj. 



For the model with t acc = tdyn, the planet is in 4: 1 resonance 
with the binary when its mass reaches and exceeds that of 
Jupiter, which arises at t - 5.5 x 10 3 P (see Fig. [2]). FigfTUl 
displays the evolution of the resonant angles associated with 
the 4:1 resonance. Comparing this figure with Fig|71 we 
see that the 4:1 resonance breaks at t ~ 8.5xl0 3 P with 
m p = 1 Mj whereas it breaks at t ~ 1.15 x 10 4 P in the one 
with m p > 1 Mj. This supports the idea that the resonant 
interaction is stronger for higher planet masses. Once again, 
the resonance drives the planet eccentricity up to e p ~ 0.2, until 
the planet undergoes a close encounter with the binary and 
is subsequently ejected from the system at t ~ 1.15 x 10 4 P. 
Here, it is worth noting that the planet mass is m p ~ 1.7 Mj 
when this occurs (see Fig|2]). 

For the model with t acc = 10 tdyn, the planet orbits just 
beyond the location of the 4:1 resonance with the binary 
when its mass attains m p = 1 Mj, which corresponds to 
t ~ 1.4 x 10 4 P. From this point in time until t ~ 2.1 x 10 4 P, 
the orbital evolution of the planet is similar to that found in 
the model with m p - 1 Mj and t acc = 10 tdyn (see Sect. 13.21) . 
Then, the large value of its eccentricity (e p ~ 0.13) causes 
the planet to undergo a close encounter with the binary and 
to be scattered out. Interestingly, this scattering leads to the 
temporary formation of the 6:1 resonance between the planet 
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Fig. 9. This figure shows the evolution of the system for models in which m p > 1 Mj. Here, the left and right panels correspond 
to t acc = tdyn and t acc = 10 td yn respectively. Top: evolution of the semi-major axis of the planet. Middle: evolution of the 
eccentricities for the planet (black), binary (red) and disc (green). Bottom: evolution of the longitudes of pericentre for the planet 
(black), binary (red) and disc (green). 



and the binary at t ~ 2.15 x 10 4 P. As can be seen in the 
middle panel of Fig{9j which shows the eccentricities of the 
planet, disc and binary, trapping into 6:1 resonance makes the 
planet eccentricity increase to e p > 0.2. At t ~ 2.3 x 10 4 P, 
a close encounter between the binary and planet occurs, and 
the planet is consequently scattered out on a high eccentricity 
orbit with a p ~ 1.8. 

Just after this scattering, the planet mass has reached 
m p ~ 3 Mj and the eccentricities of the disc and planet have 
attained ~ 0.08 and e p -0.1 respectively. In agreement with 
Nelson (2003), we find that the interaction between the eccen- 
tric disc and the eccentric planet induces outward migration of 
the latter. The total torque exerted by the disc on the planet, 
as well as the torques due to the disc interior and exterior to 
a p , are presented in Fig(T2j We see that the time-averaged total 
disc torque is clearly positive from t ~ 23 x 10 4 P onward, 
and that the torques oscillate with a period of ~ 3 x 10 3 binary 
orbits. Interestingly, these time variations correspond closely 
to the precession of the planet relative to that of the disc. 
The reason for this is simply that the eccentric disc exerts an 
orbit-averaged positive torque on the eccentric planet when 
the orbits are misaligned, and exerts a negative torque when 
the orbits are aligned. An antialigned configuration causes the 
planet to interact with matter whose angular velocity is greater 
than that of the planet at apocentre, and this leads to the planet 
experiencing a strong positive torque. 



At later times, the simulation indicates that the outward migra- 
tion of the planet can be sustained over long time scales. There 
are a number of reasons for this. First, the planet maintains 
an eccentric orbit and experiences a strong positve torque at 
apocentre when the disc and planet orbits are antialigned; this 
positive torque experienced by the planet implies a negative 
torque experienced by the outer disc material, causing gas 
to flow through the planet orbit to form an inner disc. Fig. 
ITTI shows the state of the disc, planet and binary at different 
times. The first panel shows the disc just before the planet 
enters the 4:1 resonance, and the second panel shows the 
system just after the planet has been scattered outward. The 
third and fourth panels show the increase in inner disc size 
and the outward migration of the planet. During its transition 
from outer to inner disc, the disc material continues to exert a 
positive torque on the planet, manifested through a corotation 
torque. Moreover, the existence of an inner disc provides 
another source of positive Lindblad torques and assists the 
planet in maintaining outward migration. 
Examination of the torques exerted by the disc on the planet 
(see Fig. [12]) shows that they weaken from t ~ 4 x 10 4 P. This 
is due mainly to the increase in gap size generated by the planet 
as it grows in mass, combined with the fact that the disc con- 
tains a finite reservoir of gas in our model. The increase in gap 
size clearly affects the torques due to both outer and inner discs. 



10 



A. Pierens and R.P Nelson: On the formation and migration of giant planets in circumbinary discs 



^acc ^-dyn 



t QC c=t d 



+000 6000 
Time (binary orbits) 



Fig. 7. This figure shows, for the simulation with m p = 1 Mj 
and t acc = tdyn, the evolution of the resonant angles fa, fa, fa 
and fa associated with the 4:1 resonance. 



Time (binary orbits) 



Fig. 10. This figure shows, for the simulation with m p > 1 Mj 
and t acc = tdyn, the evolution of the resonant angles fa, fa, fa 
and fa associated with the 4:1 resonance. 



At the end of the run, the planet has migrated to the outer 
edge of the disc and its mass has attained m p ~ 5 Mj. This sug- 
gests that another avenue for evolution of giant planets whose 
mass exceeds 1 Mj is long-term outward migration, in addition 
to the possibility of being scattered out of the system as dis- 
played by the runs described previously. Establishing the ratio 
of these two outcomes will require a large suite of simulations 
with slightly varying initial conditions, task which is beyond 
the scope of the paper. It appears from our results, however, 
that the probability of scattering and ejection is greater than 
that of prolonged outward migration. 

4. Summary and conclusion 

In this paper we have presented the results of hydrodynamic 
simulations aimed at studying the formation and evolution of 
giant planets embedded in circumbinary discs. 
We focused on a model in which a 20 M© core initially trapped 
at the edge of the cavity formed by the binary can slowly ac- 
crete material from the disc. We examined how the final out- 
come of the system depends on the accretion rate onto the 
planet and on the final planet mass attained. For each value of 
the accretion rate considered, we performed three calculations. 
In two of the three simulations, we assumed that accretion stops 
when the mass of the planet has reached either m p = 1 M s or 
m p = 1 Mj. In the remaining case, we allowed the planet to ac- 
crete gas freely from the disc in such a way that its final mass 
was m p > 1 Mj. The simulations show different outcomes, de- 
pending on the final mass of the planet: 



i) In models with m p = I Ms , the planet migrates inward until 
its eccentricity becomes large enough for the disc torques ex- 
erted on the planet to become positive, leading to reversal of 
migration. The planet then enters in a regime of runaway out- 
ward migration until it passes through the 5:1 resonance with 
the binary, which halts the runaway migration. From this time 
onward the planet first drifts outward very slowly, and then mi- 
grates inward again until the end of the simulation. Although 
it is not possible to run these simulations further, we speculate 
that this process of periodic inward and outward migration is 
repeated, leading to the formation of a long-term stable cir- 
cumplanetary system. 

ii) In models with final mass m p = 1 Mj, the evolution is 
found to depend weakly on the value of the accretion rate. For 
tacc = tdyn, the planet becomes locked into the 4:1 resonance 
with the binary until it is scattered to a larger radius due to 
a close encounter with the secondary star. For t acc = 10 td yn , 
the planet orbits first in, and then just outside, 4:1 resonance. 
Scattering also arises in this case. In both cases the planet mi- 
grates back in toward the binary system, and for t acc = 10 tdyn 
the planet undergoes another close encounter with the central 
binary which leads to ejection from the system. The run with 
tacc = tdyn shows slow inward migration of the planet, such that 
we are unable to run this model beyond ~ 1.6 x 10 5 binary or- 
bits. We speculate, however, that recurrent episodes of inward 
migration followed by close encounters with the central binary 
will eventually lead to ejection from the system. An alternative 
outcome, however, is scattering out to a large distance within 
the disc where the planet can orbit stably until disc dispersal. 
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Fig. 11. This figure shows, for the model in which m p > 1 Mj and t acc 
shown above the plots. Here, the planet is represented by a white circle. 



10 tdyn, snapshots of the disc surface density at times 



Once again it will require a large suite of simulations to deter- 
mine the ratios of these various outcomes, which goes beyond 
the scope of this paper. 

iii) In models with m p > 1 M/, the final outcome depends 
more strongly on the value of the accretion rate. In the sim- 
ulation with t acc = tdyn, trapping into 4:1 resonance with the 
binary leads ultimately to a planet with mass m p ~ 1.7 Mj 
being completely ejected from the system. A close encounter 
between the planet and binary is also observed in the calcula- 
tion with t acc = 10 tdyn- Subsequently however, the planet is 
not ejected but is rather scattered out to a larger radius on a 
high eccentric orbit. After this scattering, the planet migrates 
outward through its interaction with an eccentric disc, until it 
reaches the outer edge of the disc. At the end of the simulation, 
the planet mass has reached m p ~ 5 Mj. 
From an observational point of view, the results of our simu- 
lations indicate that Saturn-mass planets are probably the best 
candidates to be found in close binary systems. Higher-mass 
planets are usually found to undergo close encounters with the 
secondary star, which raises the possibility that these systems 
become unstable in the course of their evolution. 
A number of other issues remain to be resolved when consider- 
ing the early stages of planet formation in circumbinary discs. 
For example, the question of whether or not planetary cores 
can grow due to planetesimal accretion needs to be examined 



in more detail. An eccentric binary can lead to the formation 
of an eccentric disc, in such a way that planetesimal accretion 
is prevented except in the outer regions of the disc. The ability 
of planetesimals to accrete, with care being taken to simulate 
the structure of the circumbinary disc, will be the subject of a 
future paper. 

Another issue relates to the fact that we have only simulated 
a two dimensional system in this paper. The close encounters 
experienced by planets with the central binary are likely to pro- 
duce significantly inclined orbits, and the resulting disc-planet 
interaction is likely to be modified (weakened) by this, lead- 
ing to potentially different outcomes to those observed in our 
2D runs. Given the very long evolution times involved, how- 
ever, performing full 3D simulations is beyond current compu- 
tational capability. There remains scope, however, for a more 
approximate 3D treatment of this problem. 

Acknowledgements. The simulations performed in this paper were 
performed on the QMUL High Performance Computing facility pur- 
chased under the SRIF iniative. 
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